Metal-to-insulator transition and electron-hole puddle formation 
in disordered graphene nanoribbons 
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The experiment ahy observed metal-to-insulator transition in hydrogenated graphene is numeri- 
cally confirmed for actual sized graphene samples and realistic impurity concentrations. The eigen- 
states of our tight-binding model with substitutional disorder corroborate the formation of electron- 
hole-puddles with characteristic length scales comparable to the ones found in experiments. The 
puddles cause charge inhomogeneities and tend to suppress Anderson localization. Even though, 
monitoring the charge carrier quantum dynamics and performing a finite-size scaling of the local 
density of states distribution, we find strong evidence for the existence of localized states in graphene 
nanoribbons with short-range but also correlated long-range disorder. 
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The experimental observation of a disorder-induced 
metal-to-insulator transition in graphene upon hydro- 
genation [1] has triggered a vivid debate on the nature 
of this transition. For high concentrations of hydrogen 
several mechanisms of gap opening have been discussed, 
such as full sp'^ to sp^ transition, localization of sp^ ar- 
eas, or erasing of midgap states [H [sj. Graphene, on the 
other hand, is a truly two-dimensional system, and the 
one-parameter scaling theory predicts that at zero tem- 
perature any finite amount of disorder should lead to An- 
derson localization (AL) [3]. Otherwise, the existence of 
a scaling function might be questionable since the Fermi 
wavelength diverges near the charge neutrality point and 
there is no spatial scale on which the conductivity is 
much larger than e^//i (Hf. So far it seems that AL has 
not been seen in disordered graphene down to tempera- 
tures of liquid- helium [6| . This surprising result has been 
attributed to strong charge carrier density fluctuations 
that break up the sample into electron- hole puddles [ij. 
Within these puddles the local chemical potential devi- 
ates enough from the charge neutrality point to allow 
for electron or hole conductivity. Mesoscopic transport 
is then determined by activated (variable-range) hopping 
or leakage between the puddles [8]. If the formation of 
electron-hole puddles is suppressed, however, AL might 
be observed. This has been reported by quite recent ex- 
periments in double-layer graphene heterostructures [9]. 

Previous theoretical work on disordered graphene 
strongly emphasizes the difference between short- and 
long-range scatterers. While the former applies to 
the case of hydrogenation, the latter rather describes 
the effect of charged impurities in the substrate [Toj . 
Within the Dirac approximation, only short-range im- 
purities cause intervalley scattering, and thus may lead 
to AL [ll| . The presence of long-range impurities alone 
gives rise to intravalley scattering which is not sufficient 
to localize the charge carriers [12]. Another factor is the 
edge geometry of the graphene nanoribbons (GNRs) that 
determines the universality class of disordered samples 




FIG. 1. (Color online) Cartoon of the substitutional disorder 
model describing hydrogenated graphene. 



as long as the phase coherence length exceeds the system 
size 



13j . Going beyond the Dirac approximation and de- 



scribing graphene by a tight-binding model, it is natural 
to ask whether the scattering range is still decisive. Devi- 
ations from the idealized linear dispersion, a finite lattice 
spacing, and the trigonal lattice symmetry, which breaks 
the rotational symmetry of the Dirac cones, call for a 
careful numerical analysis of the localization properties 
within the tight-binding description [14| . 

In this work, we prove by unbiased numerics that ex- 
perimentally relevant concentrations of hydrogen x <C 1 
may induce a metal-to-insulator transition in actual- 
size graphene samples. We show that the single-particle 
wavefunctions of our disorder model are also localized for 
correlated long-range disorder. Even for potential fluc- 
tuations on an atomistic scale there is strong evidence 
for electron-hole puddle formation on an intrinsic scale 
of some ten nm, in agreement with recent experimen- 
tal observations. [l5|. In contrast to previous studies on 
temperature dependent transport in disordered graphene 
within the semiclassical Boltzmann approach [16], we re- 
strict ourselves in the following to strictly zero tempera- 
ture and adopt a purely quantum point of view. 

We consider a tight-binding Hamiltonian H = 
-t Y^^ij) {c\cj + H.c.) + Y.. Vic\c^, on the honeycomb lat- 
tice with TV sites, where the operators c\ (q) create 
(annihilate) an electron in a Wannier state centered at 
site and t denotes the nearest-neighbor transfer in- 
tegral. The landscape of onsite potentials {Vi} results 
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FIG. 2. (Color online) Spatial distribution of the normalized 
LDOS pi/ pme E — for the binary alloy model with poten- 
tial difference A = 6.0t and impurity concentration x — 0.1% 
(left column) and 1% (right column). The disorder configura- 
tions shown in the top panel were used for both BC. Data 
obtained by exact diagonalization (ED), GNR sample size 
(37 X 12) nm^, corresponding to 300x60 atoms. 



from the superposition of contributions of A^imp = 
randomly distributed Gaussian impurities at positions 
V, = E™=rA^exp(-|ri-r™|V(2^2)) gy 

choice of ^ the range of the individual impurity potentials 
can be continuously tuned from short-ranged to long- 
ranged. For ^ ^ and x = 1 we recover the Ander- 
son model on a GNR [l8[. Assuming a fixed = A 
for all impurities, the limit ^ — > results in the bi- 
nary alloy model in which only distinct sites acquire 
a finite onsite potential. Vacancies correspond to sites 
with A ^ oo, leading to a quantum site-percolation sce- 
nario [HI. The presence of adsorbed hydrogen atoms 
alters the hybridization of carbon atoms from sp^ to sp^, 
partially removing the corresponding Pz orbital from the 
TT-band. We model the yet finite probability of finding 
electrons at the adsorbant site by a finite value of the 
disorder strength A (see Fig. [T|). 

Experimental results by Bostwick et al. [H suggest a 
metal-to-insulator transition in graphene for a hydrogen 
coverage as low as 0.3%. In Fig. [2] we contrast the spa- 
tial distribution of the local density of states (LDOS) 
pi{E) = \{n\i)\'^6{E — En) at the Dirac point energy 
for a hydrogen coverage slightly below and above this 
threshold. For zigzag boundaries the well known edge 
states persist even in the presence of weak disorder. Im- 
purities on the sublattice with high LDOS values near 
the GNR edges drastically reduce these values. On the 
other sublattice they do not have any effect. In the bulk 
of the ribbon, the LDOS is slightly enhanced as com- 
pared to the ordered case. Positive interference traps the 
wavefunction on sites in between the impurities. For pe- 
riodic boundary conditions (PBC) the spatial LDOS dis- 
tribution is clearly distinct for both impurity concentra- 
tions: We observe only slight local perturbations of the 
perfectly extended state for low impurity concentrations 
but a clearly localized state at 1% hydrogen coverage. 
Measurements on a sample of this size therefore yield 
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FIG. 3. (Golor online) Upper panel: Distribution of the LDOS 
at experimentally relevant energies for the binary alloy model 
with different impurity concentrations and PBG. The sample 
width W = 109 nm. Normalization of the distribution to pme 
directly relates its position to the height of the maximum and 
its width. Inset: Magnification of the averaged DOS with 
indications of the energies for which the LDOS distributions 
are shown. Lower panel: Finite-size scaling of the LDOS 
distribution. Data obtained by the kernel polynomial method 
with resolution adapted to the level-spacing, Nk = 140 (for 
details see fU,!!!). 



metallic (insulating) behavior for coverages 0.1% (1%). 
The observed metallic character seems to be merely a 
consequence of a finite localization length A that exceeds 
the system size. Note that the calculated state charac- 
teristics and experimental results agree qualitatively for 
PBC only. This underlines that the observed localization 
properties are intrinsic to short-range disordered bulk 
graphene. Edge effects arise on top, but are to a certain 
extent irrelevant in experiments, especially if mobilities 
are measured using a multiterminal Hall geometry (2q| . 

In order to assert that AL takes place in infinite GNRs, 
we analyze the distribution of the normalized LDOS in 
Fig. [3l We restrict ourselves to three characteristic en- 
ergies which are shown in the inset together with the 
averaged density of states pme = (pi)- The behavior 
of the LDOS distribution upon finite-size scaling (lower 
panel) is a powerful criterion to detect AL for different 
kinds of disordered systems even in presence of in- 
teractions [23]. Extended states are characterized by an 
f [pi /pme] being independent of the system size. Other- 
wise sensitivity of the distribution to the system size in- 
dicates localization, which we indeed observe for binary 
alloy disordered GNRs for all energies and both impurity 
concentrations. 

For a given state, the shape of the LDOS distribution 
and the extent of its shifting depend on A; the more pro- 
nounced the shift and the more asymmetric f [pi /pme], 
the shorter is A. Larger impurity concentrations enhance 
localization, as can be seen from the asymmetric shape 
of f[pi/pme] for X = 0.01. The persisting size depen- 
dence for X = 0.001 proves localization also for such a 
weak randomness. Here the high sensibility of the LDOS 
distribution to the ratio of A and system size is of vital 
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FIG. 4. (Color online) Time evolution of the local particle 
density on disordered zigzag GNRs after, at r = 0, finite (or- 
dered) leads were attached to the left and right of the GNR. 
|'0(ri,t)|^ is normalized to the actual mean particle density 
on the GNR. The color scale is identical to Fig. [2l Sample 
dimensions are (221 x 109) nm^, corresponding to 1800 x 512 
lattice sites. Left and middle column: Binary alloy model 
with A = 6t. Right column: Gaussian correlated disorder 
model with ^ = 3a, where a is the inter-carbon distance. 
Here the potential is normalized to max(Vi) = 6t. The in- 
sets show magnifications of \ip(ri,t)\'^ in the quasistationary 
regime together with the corresponding potential landscape. 



x = 0.001 





FIG. 5. (Color online) Quasistationary local particle density 
summed over the transverse direction. Each subpanel shows 
how changing one control parameter influences the baseline 
case (red dashed line: binary disorder, x — 0.01, E — 0, 
zigzag BC) while keeping the others fixed. 



importance. It allows to detect localization also in the 
case of weak disorder for which A distinctly exceeds the 
system size and consequently /[pi/pme] is concentrated 
around unity. 

In Fig. m we contrast the quantum dynamics of a par- 
ticle injected into zigzag GNRs with binary or Gaussian 
correlated disorder and impurity concentrations of 0.1% 
and 1%. As initial state |V^o) we prepare an exact E = {) 
eigenstate of the ordered infinite graphene lattice in the 
lead left to the sample (see Fig. [T]). After bringing the 
lead in contact with the sample, we let the system evolve 
in time by solving the time-dependent Schrodinger equa- 
tion by a Chebyshev expansion technique [TqI, . Due 
to the coupling with the disordered GNR, |?/^o) is not 
an eigenstate of the overall system but comprises admix- 
tures of other states, mainly from the vicinity of £^ = 0. 
The snapshot at r = IO^tq, where tq = h/t^ confirms the 



intuition that spreading is faster the lower the impurity 
concentration is (see Fig. [4j). At r = IO^tq all states 
have reached quasistationarity, and we can extract their 
characteristic features. For x = 0.1% the state spans the 
whole sample. The inset reveals its puddle-like structure 
with density fluctuations of two orders of magnitude on 
length scales of 5 — 10 nm. At an impurity concentra- 
tion of X = 1% the particle density is reduced about two 
orders of magnitude between left and right edge of the 
GNR, providing a direct visualization of AL. Note that 
the local structure of |?/^) remains puddle-like, but the 
spatial extent of the puddles is substantially reduced be- 
low 1 nm. From the similarity of both local structures one 
may argue that for larger systems also states for x = 0.1% 
will be localized. Note that correlated disorder results in 
a markedly smoother potential landscape. In addition, 
the electron-hole puddles are superimposed by a coarse- 
grained filamentary structure. 

Integrating the local particle density over the trans- 
verse ribbon direction allows for a more quantitative 
analysis of the localization properties (see Fig. [SJ. For a 
fixed configuration of impurity positions we vary the rel- 
evant control parameters and extract A from fitting the 
quasistationary density to an exponential decay. Remov- 
ing part of the impurities results in a larger A [Fig.[5ja)]. 
For the considered ribbon width the difference between 
zigzag- and PEG is marginal, while armchair edges dras- 
tically reduce the transmission [Fig. [5jb)]. This reflects 
the mismatch of preferred transport direction and ribbon 
axis. Varying the incident particle energy [Fig. Efc)] we 
observe an enhanced transmission for E = — 0.2t, which 
might be attributed to resonant (localized) states. Note 
that at the position of the chemical potential of graphene 
on a SiC substrate, E = 0.2t, A is even shorter than at 
the Dirac point ^ = 0. For correlated long-range dis- 
order [Fig. Efd)], A increases with ^, yielding extended 
states if ^ a. For ^ = a, 3a the states are still clearly 
localized. This disagrees with results in the literature 
obtained within the Dirac approximation [l2| , where the 
valley quantum number is conserved and localization is 
suppressed in absence of intervalley scattering ^Jj] . We 
attribute this difference to the lattice discreteness and the 
breaking of the rotational symmetry of the Dirac cones by 
the trigonal symmetry of the honeycomb lattice. More- 
over, even for a narrow banded initial state the dynamics 
is influenced by states from the whole energy spectrum 
where the graphene dispersion significantly deviates from 
the linear approximation. If all these aspects were taken 
into account, long-ranged disorder may cause localiza- 
tion within a tight-binding description in accordance with 
one-parameter scaling [25| . Whether AL really occurs for 
correlated disorder can be proven by performing a finite- 
size scaling of the LDOS distribution. In doing so we 
find evidence for localization for ^ = 3a from the shift- 
ing of the LDOS distribution (see left panel of Fig. [6|). 
A remarkable feature of /[pi/pme] is the salient tail that 
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FIG. 6. (Color online) Left: Finite-size scaling of the LDOS 
distributions for the Gaussian correlated disorder model with 
parameters matching Fig. [H Right: CDF of the LDOS for 
various potential ranges ^. 




FIG. 7. (Color online) Spatial distribution of the normalized 
LDOS at ^ = for the disorder configuration shown in the 
left panel with PBC. System sizes and color mapping are the 
same as in Fig. O Data obtained by ED. 

develops for small values of pi on increasing ^. In the 
right panel of Fig. [6] this results in a kinked cumulated 
distribution function (CDF) instead of the approximate 
straight CDF for uncorrelated disorder (26| . 

So far we considered strong disorder (A ^ 1) and 
weak randomness {x <C 1), for which it is tempting to 
relate average puddle size and distance between the im- 
purities. Interestingly, electron- hole puddles also arise 
for weak disorder strength and strong randomness, corre- 
sponding to a disorder landscape varying on an atomistic 
scale without any correlations. Such a modeling might 
be regarded as an attempt to capture the effect of the 
buffer layer forming between epitaxially grown graphene 
and its SiC substrate \^\. In this setup we refrain from 
considering a correlated potential landscape in order not 
to impose any a priori correlations in the LDOS. For weak 
disorder strength, A = 0.5t, the LDOS at £^ = never- 
theless becomes puddle-like with a characteristic scale of 
2-5 nm (see middle panel of Fig. [7j). Note that the choice 
of BC has no qualitative impact on the LDOS in this 
limit of the binary alloy model. Obviously, any subtle 
differences in the state characteristics induced by BC are 
masked by the randomness of the potential landscape. 
Increasing the potential difference to A = 2t, A drops 
below the system size and we observe clearly localized 
states. 

To conclude, hydrogenated graphene behaves at zero 
temperature as a 'normal' two-dimensional disordered 
system concerning AL, provided the extensions of ultra- 
high-quality samples become very large. If the localiza- 
tion length noticeably exceeds the system size, the sample 
nevertheless shows metallic behavior. We find that also 
certain long-range correlations in the disorder landscape 
yield localized single-particle wavefunctions. Most no- 
tably, we show that disorder-induced electron-hole pud- 



dles may arise for both disorder types. The intrinsic scale 
of the puddle-like structures in the eigenstates is not 
simply set by the distance between impurities, but re- 
sults from subtle quantum interference effects. Even for 
atomic scale fluctuations of the disorder potential they 
might exceed 5 nm, which is in the range of experimen- 
tally measured values. The presence of electron- hole pud- 
dles, leading to intra- and inter-puddle transport, drives 
the system away from the metal-to-insulator transition, 
thereby masking AL [9] . This resolves the puzzle why AL 
is so hard to detect in disordered graphene and GNRs. 
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